// QSFLAME.C - "cosmic recursive fractal flames"
// GNU version by Scott Draves <http://scottdraves.com> Copyright (C) 1992
// W95 version by Michael Sargent <msargent@zoo.uvm.edu> Copyright (C) 1999
// Revisions by Michael Sargent <msargent@uvm.edu> Copyright (C) 2018
//
// This file contains the flame-specific functions and WinMain()
//
// Needed for compiling: this file
//                       QSFGUI.C
//                       QSFLAME.H
//                       QSFLAME.RC
//                       SPLASH1.BMP  - 320 x 320
//                       SPLASH2.BMP  - 320 x 320
//                       ABOUT.BMP    - 480 x 240
//                       QSFLAME.ICO  -  32 x  32
/*=========================================================================
   This is unsupported, non-standard software. It is provided "as is" and
   any express or implied warranties, including, but not limited to, the
   implied warranties of merchantability and fitness for a particular
   purpose are disclaimed. In no event shall the author be liable for any
   direct, indirect, incidental, special, exemplary, or consequential
   damages (including, but not limited to, procurement of substitute goods
   or services; loss of use, data or profits; or business interruption)
   however caused and on any theory of liablility, whether in contract,
   strict liability, or tort (including negligence or otherwise) arising in
   any way out of the use of this software, even if advised of the
   possibility of such damage.
=========================================================================
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, see <https://www.gnu.org/licenses/>.
//=======================================================================*/
#include "QSFlame.h"
#include "QSFGUI.c"
//=========================================================================
// WinMain
//=========================================================================
#pragma argsused
int WINAPI WinMain(HINSTANCE hInstance, HINSTANCE hPrevInst,
                   LPSTR lpszArgs, int nCmdShow)
{
   int error;

   //----- initialize and register the window classes -----
   error = InitWindowClass(hInstance);

   //----- create the main window -----
   if (!error) error = CreateMainWindow(hInstance, nCmdShow);

   //----- call the major program functions and run the message loop -----
   if (!error) return ProgramDriver();
   else
   {
      MessageBox(NULL, errormessage[error - 1], "Error",
                 MB_OK | MB_ICONEXCLAMATION);
      return 0;
   }
}
//=========================================================================
// Set up the color map array
// At this point, p.white_level is held constant at 255, so this function
//    only copies the palette[][] array to the cmap[][] array and sets
//    all cmap[j][3] values to 255
//=========================================================================
void setup_cmap(void)
{
   int j, k;
   double cfactor;

   //----- compute the colormap entries -----
   cfactor = (double)p.white_level / 255.0;
   for (j = 0; j < 256; j++)
   {
         for (k = 0; k < 3; k++)
            cmap[j][k] = (int)((double)palette[j][k] * cfactor);
         cmap[j][3] = p.white_level;
   }
}
//=========================================================================
// Set up a simplified Gaussian filter; this is a streamlined version of
//    Scott Draves's original formula - the filter is stored in the linear
//    integral table "itable[]" and is applied first vertically and then
//    horizontally by the apply_filter() function
// Should only be called in truecolor modes and if fbufferflag == TRUE
// "radius" is the filter radius in pixels; thus it is an integer
//=========================================================================
void setup_filter(int radius)
{
   int i;
   double ii, jj, sum;

   jj = ((2.0 * (double)radius + 1.0) /
         (double)((radius << 1) + 1) - 1.0) * FILTER_CUTOFF;
   for (i = 0; i < (radius << 1) + 1; i++)
   {
      ii = ((2.0 * (double)i + 1.0) /
            (double)((radius << 1) + 1) - 1.0) * FILTER_CUTOFF;
      itable[i] = exp(-2.0 * (ii * ii + jj * jj));
      if (i > 0) itable[i] += itable[i - 1];
   }
   sum = itable[i - 1];
   for (i = 0; i < (radius << 1) + 1; i++) itable[i] /= sum;
}
//=========================================================================
// Set up a gamma lookup table using the standard gamma correction formula
// Should only be called in truecolor modes
//=========================================================================
void setup_gamma(double gamma)
{
   int i;
   double rgamma;

   rgamma = 1.0 / gamma;
   for (i = 0; i < 256; i++)
      gtable[i] = (unsigned char)(256.0 * pow((double)i / 256.0, rgamma));
}
//=========================================================================
// Calculate and plot (p.iter == 100,000) points
// Set up prior to the iteration loop:
//    clear the display and virtual screens
//    set up the "camera" dimensions in abstract units
//    set up the intensity factor calculations
//    the initial point is randomized in "point[]", or if this is a
//       "Re-iteration" the previous loop's last point is picked up.
//    set up the xform_distrib[] array
//    zero the buffers if this is a new image
//=========================================================================
void iterate(void)
{
   HDC hdc;                                 // Windows GUI variables
   RECT CurrentClientWindow;
   HBRUSH hBrush;
   HCURSOR hArrow, hHourglass;
   MSG msg;

   COLORREF cr;                             // holds GetPixel() return value
   int i, j, k, startindex;                 // loop counters
   int offset;                              // for reading fbuffer[]
   int temp;                                // for setting brightness values
   int xform_distrib[CHOOSE_XFORM_GRAIN];   // array which holds the indices
                                            //    of the active xforms, in
                                            //    the ratio of their "density"
                                            //    factors
   int fn;                                  // randomly-chosen index for
                                            //    xform_distrib[]
   int h, w;                                // x/y location of pixels in the
                                            //    virtual image, which is
                                            //    larger than the final image
                                            //    if "oversample" is > 1
   int color_index;                         // derived from point[2] as the
                                            //    initial index for cmap[][]
                                            //    for each calculated pixel
   int channel[4];                          // stores revised RGBA values
                                            //    for each pixel
   double c[4];                             // floating point mirror of the
                                            //    channel[] array
	double ppu, sample_density;              // values from "p" struct,
                                            //    adjusted by p.scale
   double area, k1, k2, ls;                 // for calculating pixel intensity
   double size[2];                          // dimensions of image in abstract
                                            //    units; used to convert a
                                            //    point's location to pixels
   double point[3];                         // location and color value of
                                            //    each calculated point in
                                            //    abstract units
   static double tpoint[3];                 // saves last point values for the
                                            //    "Re-Iterate" option
   double t, r, dr;                         // for setting up xform_distrib[]
   double tx, ty, v, nx, ny, r2, c1, c2, a; // iteration variables
 
	hdc = GetDC(g_hWnd);

   //----- clear the display and virtual screens -----
   if (clrscr)
   {
      hBrush = (HBRUSH)GetStockObject(BLACK_BRUSH);
      SelectObject(hdc, hBrush);
      GetClientRect(g_hWnd, &CurrentClientWindow);
      FillRect(hdc, &CurrentClientWindow, hBrush);
		hBrush = CreateSolidBrush(RGB(rback, gback, bback));
      SelectObject(hdc, hBrush);
      CurrentClientWindow.right = min(xres, ClientRes.right);
      CurrentClientWindow.bottom = min(yres, ClientRes.bottom);
      FillRect(hdc, &CurrentClientWindow, hBrush);
      SelectObject(g_vhdc, hBrush);
      CurrentClientWindow.right = xres;
      CurrentClientWindow.bottom = yres;
      FillRect(g_vhdc, &CurrentClientWindow, hBrush);
      DeleteObject(hBrush);
   }

   //----- set up the "camera" dimensions in abstract units -----
   ppu = p.pixels_per_unit * p.scale;
   sample_density = p.sample_density * p.scale * p.scale;
   bounds.left = p.center[0] - 4.0 / p.scale;
   bounds.top = p.center[1] - 3.0 / p.scale;
   bounds.right = bounds.left + 8.0 / p.scale;
   bounds.bottom = bounds.top + 6.0 / p.scale;
   size[0] = 1.0 / (bounds.right - bounds.left);
   size[1] = 1.0 / (bounds.bottom - bounds.top);

   //----- set up the intensity factor calculations -----
   area = ((double)xres * (double)yres) / (double)(ppu * ppu);
   k1 = (p.contrast * (double)PREFILTER_WHITE * 268.0) / 256.0;
   k2 = (double)(oversample * oversample) /
   	(p.contrast * area * (double)p.white_level * sample_density);

   //----- if there is no hit buffer, "ls" only needs to be -----
   //-----    calculated once -----
   if (hitflag == FALSE)
   	ls = (k1 * log(1.0 + (double)p.white_level * k2)) / (double)p.white_level;

   //----- initial point -----
   //----- pick up previous loop's last point if reiterating -----
   if (reitflag)
      for (j = 0; j < 3; j++) point[j] = tpoint[j];
   else
   {
      point[0] = qsrandom11();
      point[1] = qsrandom11();
      point[2] = qsrandom01();
   }

   //----- "set up xform_distrib[], which is an array that converts -----
   //-----    a uniform random variable into one with the distribution -----
   //-----    dictated by the density fields" (S.D.) -----
   dr = 0.0;
   for (j = 0; j < NXFORMS; j++) dr += x[j].density;
   dr = dr / CHOOSE_XFORM_GRAIN;
   j = 0;
   t = x[0].density;
   r = 0.0;
   for (k = 0; k < CHOOSE_XFORM_GRAIN; k++)
   {
      while (r >= t)
      {
	      j++;
	      t += x[j].density;
      }
      xform_distrib[k] = j;
      r += dr;
   }

   //----- zero the buffers if this is a new image -----
   if (hitflag && clrscr) ZeroMemory(hits, xos * yos * sizeof(short));
   if (fbufferflag && clrscr)
   {
      for (i = 0; i < xos * yos * 3; i += 3)
      {
         fbuffer[i] = (unsigned char)rback;
         fbuffer[i + 1] = (unsigned char)gback;
         fbuffer[i + 2] = (unsigned char)bback;
      }
   }

   //----- display hourglass mouse cursor while iterating -----
   hArrow = LoadCursor(NULL, IDC_ARROW);
   hHourglass = LoadCursor(NULL, IDC_WAIT);
   SetCursor(hHourglass);

   //----- iteration loop; p.iter == 100,000 -----
   //----- if this is a new image, ignore the first FUSE iterations -----
   if (reitflag) startindex = 0;
   else startindex = -FUSE;
   for (i = startindex; i < p.iter; i++)
   {
      //----- check for explosion -----
      if ((fabs(point[0]) > 1E6) || (fabs(point[1]) > 1E6)) break;

      fn = xform_distrib[random(CHOOSE_XFORM_GRAIN)];

      //----- "first compute the color coord..." -----
      point[2] = (point[2] + x[fn].color) / 2.0;

      //----- "... then apply the affine part of the function..." -----
      tx = x[fn].c[0][0] * point[0] + x[fn].c[1][0] * point[1] +
	        x[fn].c[2][0];
      ty = x[fn].c[0][1] * point[0] + x[fn].c[1][1] * point[1] +
	        x[fn].c[2][1];
      point[0] = point[1] = 0.0;

      //----- "... then add proportional amounts of each -----
      //-----    variation" (S.D.) -----

      //----- LINEAR -----
	   v = x[fn].var[0];
	   if (v > 0.0)
	   {
	      nx = tx;
	      ny = ty;
	      point[0] += v * nx;
	      point[1] += v * ny;
	   }

      //----- SINUSOIDAL -----
	   v = x[fn].var[1];
	   if (v > 0.0)
	   {
	      nx = sin(tx);
	      ny = sin(ty);
	      point[0] += v * nx;
	      point[1] += v * ny;
	   }

      //----- SPHERICAL -----
	   v = x[fn].var[2];
	   if (v > 0.0)
	   {
	      r2 = tx * tx + ty * ty + 1e-6;
	      nx = tx / r2;
	      ny = ty / r2;
	      point[0] += v * nx;
	      point[1] += v * ny;
	   }

      //----- SWIRL -----
	   v = x[fn].var[3];
	   if (v > 0.0)
	   {
	      r2 = tx * tx + ty * ty;    // "k here is fun" (S.D.)
	      c1 = sin(r2);
	      c2 = cos(r2);
	      nx = c1 * tx - c2 * ty;
	      ny = c2 * tx + c1 * ty;
	      point[0] += v * nx;
	      point[1] += v * ny;
	   }

      //----- HORSESHOE -----
	   v = x[fn].var[4];
	   if (v > 0.0)
	   {
	      if (tx < -EPS || tx > EPS || ty < -EPS || ty > EPS)
	         a = atan2(tx, ty);      // "times k here is fun" (S.D.)
	      else
	         a = 0.0;
	      c1 = sin(a);
	      c2 = cos(a);
	      nx = c1 * tx - c2 * ty;
	      ny = c2 * tx + c1 * ty;
	      point[0] += v * nx;
	      point[1] += v * ny;
	   }

      //----- POLAR -----
	   v = x[fn].var[5];
	   if (v > 0.0)
	   {
	      if (tx < -EPS || tx > EPS || ty < -EPS || ty > EPS)
	         nx = atan2(tx, ty) / M_PI;
	      else
	         nx = 0.0;
	      ny = sqrt(tx * tx + ty * ty) - 1.0;
	      point[0] += v * nx;
	      point[1] += v * ny;
	   }

      //----- BENT -----
	   v = x[fn].var[6];
	   if (v > 0.0)
	   {
	      nx = tx;
	      ny = ty;
	      if (nx < 0.0) nx = nx * 2.0;
	      if (ny < 0.0) ny = ny / 2.0;
	      point[0] += v * nx;
	      point[1] += v * ny;
	   }

      //----- if over FUSE, plot the point -----
      if (i >= 0)
      {
         if (hitflag) channel[3] = p.white_level;

         //----- if the point is within the bounds of the image, -----
         //-----    get its x/y coordinates in the virtual image -----
         //-----    and get its color_index value -----
	      if (point[0] < bounds.left || point[1] < bounds.top ||
             point[0] > bounds.right || point[1] > bounds.bottom) continue;
	      w = (int)(xos * (point[0] - bounds.left) * size[0]);
	      h = (int)(yos * (point[1] - bounds.top) * size[1]);
         if (w < 0) w = 0;
         else if (w > xos - 1) w = xos - 1;
         if (h < 0) h = 0;
         else if (h > yos - 1) h = yos - 1;
	      color_index = (int)(point[2] * 256.0);
         if (color_index < 0) color_index = 0;
	      else if (color_index > 255) color_index = 255;

         //----- adjust and plot the 24-bit pixel, and store -----
         //-----    it in the image buffer for filtering -----
         //----- This option is my attempt to get as close as I -----
         //-----    could to Scott's coloring algorithm while -----
         //-----    minimizing humongous buffers. It uses one -----
         //-----    accumulator buffer and an additional "hit -----
         //-----    buffer" of (xos * yos * sizeof(short int)) -----
         //-----    bytes to store the number of times any pixel -----
         //-----    is hit. While we are accumulating processed -----
         //-----    rather than raw pixels, the results from -----
         //-----    both approaches turn out to be so similar that -----
         //-----    it is hard to justify the additional code and -----
         //-----    processing time required by the more -----
         //-----    meticulous approach. -----
         //----- channel[] holds the working values of the RGBA ----
         //-----    attributes of each pixel as they are adjusted; -----
         //-----    it is type int so it can temporarily exceed 255 -----
         if (hitflag)
         {
         	//----- easier to read buffer than virtual screen -----
            if (fbufferflag)
            {
            	offset = (linetable[h] + w) * 3;
               channel[0] = (int)fbuffer[offset];
               channel[1] = (int)fbuffer[offset + 1];
               channel[2] = (int)fbuffer[offset + 2];
            }
            else
            {
            	cr = GetPixel(g_vhdc, h, w);
               channel[0] = (int)GetRValue(cr);
               channel[1] = (int)GetGValue(cr);
               channel[2] = (int)GetBValue(cr);
            }

            //----- get the cumulative value of the pixel so far in -----
            //-----    the loop, add the current value, and increment -----
            //-----    the hit buffer -----
            for (k = 0; k < 4; k++)
            {
            	channel[k] *= hits[linetable[h] + w];
               channel[k] += cmap[color_index][k];
               c[k] = (double)channel[k];
            }
            ++hits[linetable[h] + w];

            //----- "log intensity in hsv space" (S.D.) -----
            //----- then convert from double to integer and divide -----
            //-----    the cumulative value by the hit count to get -----
            //-----    channel[] values between 0 and 255 -----
            if (c[3] <= 0.0) continue;
            ls = (k1 * log(1.0 + c[3] * k2)) / c[3];
            for (k = 0; k < 4; k++)
            {
            	c[k] += (c[k] * ls) + 0.5;
               channel[k] = (int)(c[k] / (double)(hits[linetable[h] + w] + 0.5));
               if (channel[k] > 255) channel[k] = 255;
            }
         }

         //----- if the hit buffer could not be allocated, this -----
         //-----    is a simpler/faster and acceptable approach -----
         else
         {
         	for (k = 0; k < 4; k++) channel[k] = cmap[color_index][k];

            //----- "log intensity in hsv space" (S.D.) -----
            for (k = 0; k < 4; k++)
            {
            	c[k] *= ls;
               bump_no_overflow(channel[k], (int)(c[k] + 0.5));
               if (channel[k] > 255) channel[k] = 255;
            }
         }

         //----- apply brightness, gamma and "alpha"  -----
         //----- plot the image to the virtual screen -----
         //----- store it in the filter buffer -----
         for (k = 0; k < 3; k++)
         {
         	if (p.brightness != 0)
            {
            	temp = channel[k] + p.brightness;
               if (temp < 0) temp = 0;
               else if (temp > 255) temp = 255;
               channel [k] = temp;
            }
            if (p.gamma != 1.0) channel[k] = gtable[channel[k]];
            channel[k] += ((channel[k] * (256 - channel[3])) >> 8);
            if (channel[k] > 255) channel[k] = 255;
         }

         //----- finally plot the pixel to the virtual screen and -----
         //-----    store it in the filter buffer -----
         SetPixel(g_vhdc, w / oversample, h / oversample,
                     RGB((BYTE)channel[0], (BYTE)channel[1],
                     (BYTE)channel[2]));
         if (fbufferflag)
         {
         	offset = (linetable[h] + w) * 3;
            for (k = 0; k < 3; k++) fbuffer[offset + k] = (BYTE)channel[k];
         }

         //----- copy image to display screen every 8092 cycles -----
         if (((i & 0x1fff) == 0) || (i == p.iter - 1))
	         BitBlt(hdc, 0, 0, xscreen, yscreen, g_vhdc, 0 + xCurrentScroll,
                      0 + yCurrentScroll, SRCCOPY);
      }

      //----- poll for keypress or mouse click to break out of loop -----
      if (PeekMessage(&msg, g_hWnd, WM_KEYFIRST, WM_MOUSELAST, PM_NOREMOVE))
      {
         TranslateMessage(&msg);
         PeekMessage(&msg, g_hWnd, WM_KEYFIRST, WM_MOUSELAST, PM_REMOVE);
         if (msg.message == WM_KEYDOWN || msg.message == WM_SYSKEYDOWN ||
         	 msg.message == WM_LBUTTONDOWN || msg.message == WM_RBUTTONDOWN)
         {
            if (infinite) infinite = FALSE;
            break;
         }
      }
   }
   for (j = 0; j < 3; j++) tpoint[j] = point[j];
	ReleaseDC(g_hWnd, hdc);
   SetCursor(hArrow);
}
//=========================================================================
// Apply modified Gaussian filter to the image
// Should only be called in truecolor modes and if fbufferflag == TRUE
//=========================================================================
void apply_filter(void)
{
   int i, j, k, m, n, count, iflag, offset, iover, jover;
   double previous, startarea, endarea;
   RGBPIXEL color, integral, temp;
   RGBPIXEL *colintegral;              // array with "xos" members
   COLORREF cr;
   HDC hdc;
   MSG msg;
   HCURSOR hArrow, hHourglass;

   //----- buffer is too small for this to fail  :-)  -----
   colintegral = (RGBPIXEL *)LocalAlloc(LPTR, xos * sizeof(RGBPIXEL));

   hdc = GetDC(g_hWnd);
   hArrow = LoadCursor(NULL, IDC_ARROW);
   hHourglass = LoadCursor(NULL, IDC_WAIT);
   SetCursor(hHourglass);

   for (j = 0; j < yres; j++)
   {
      jover = j * oversample;

      //----- this loop integrates the filter for each pixel in the -----
      //-----    current horizontal line (j) along a vertical column -----
      //-----    of size (2 * radius + 1) -----
      //----- if oversample > 1, this has to be done for each -----
      //-----    subpixel (hence the nested "k" loop) -----
      for (i = 0; i < xres; i++)
      {
         for (k = 0;  k < oversample; k++)
         {
            offset = xres * k;
            colintegral[i + offset].r = colintegral[i + offset].g =
               colintegral[i + offset].b = 0.0;

            //----- integrate in each column over (2 * radius + 1) -----
            //-----    pixels, in the legal range of j, and store -----
            //-----    results in the colintegral[] array -----
            for (n = jover - p.radius, count = iflag = 0,
                 previous = 0.0; n < jover + p.radius + 1; n++)
            {
               if (n >= 0 && n < yos)
               {
                  getRGBPIXEL(&color, i + offset, n);
                  if (!iflag)
                  {
                     iflag = 1;
                     startarea = previous;
                  }
                  RGBMULT(color, (itable[count] - previous));
                  RGBADD(colintegral[i + offset], color);
               }
               previous = itable[count];
               if (n < yos) count++;
            }
            endarea = 1.0 / (itable[count - 1] - startarea);
            RGBMULT(colintegral[i + offset], endarea);
            RGBCLAMP(colintegral[i + offset]);
         }
      }

      //----- now this loop integrates the filter for each pixel in the -----
      //-----    current horizontal line along a horizontal sub-row -----
      //-----    of size (2 * radius + 1) -----
      //----- if oversample > 1, "colintegral[]" has  been calculated -----
      //-----    in subpixels by the preceding loop -----
      for (i = 0; i < xres; i++)
      {
         iover = i * oversample;
         integral.r = integral.g = integral.b = 0.0;

         //----- add contribution of each column integral in the -----
         //-----    legal range of i -----
         for (m = iover - p.radius, count = iflag = 0, previous = 0.0;
              m < iover + p.radius + 1; m++)
         {
            if (m >= 0 && m < xos)
            {
               temp = colintegral[m];
               if (!iflag)
               {
                  iflag = 1;
                  startarea = previous;
               }
               RGBMULT(temp, (itable[count] - previous));
               RGBADD(integral, temp);
            }
            previous = itable[count];
            if (m < xos) count++;
         }
         endarea = 1.0 / (itable[count - 1] - startarea);
         RGBMULT(integral, endarea);
         RGBCLAMP(integral);

         cr = RGB((unsigned char)(integral.r),
                  (unsigned char)(integral.g),
                  (unsigned char)(integral.b));
         SetPixel(g_vhdc, i, j, cr);
      }

      //----- bit-blit to the display line-by-line -----
      if (j >= yCurrentScroll && j < lastdisplayline)
         BitBlt(hdc, 0, j - yCurrentScroll, ClientRes.right, 1,
                g_vhdc, 0 + xCurrentScroll, j, SRCCOPY);

      //----- poll for keypress to break out of loop -----
      if (PeekMessage(&msg, g_hWnd, WM_KEYFIRST, WM_KEYLAST, PM_NOREMOVE))
      {
         TranslateMessage(&msg);
         PeekMessage(&msg, g_hWnd, WM_KEYFIRST, WM_KEYLAST, PM_REMOVE);
         if (msg.message == WM_KEYDOWN || msg.message == WM_SYSKEYDOWN) break;
      }
   }
   ReleaseDC(g_hWnd, hdc);
   LocalFree((HLOCAL)colintegral);

   //----- if oversample == 1, copy filtered image to fbuffer; -----
   //-----    otherwise, the HDC and the buffer don't correspond -----
   //----- didn't use GetDIBits() because it reverses the RGB order -----
   if (oversample == 1)
   {
      count = 0;
      for (j = 0; j < yres; j++)
         for (i = 0; i < xres; i++)
         {
            cr = GetPixel(g_vhdc, i, j);
            fbuffer[count++] = GetRValue(cr);
            fbuffer[count++] = GetGValue(cr);
            fbuffer[count++] = GetBValue(cr);
         }
   }
   SetCursor(hArrow);
}
//=========================================================================
// Return uncorrected pixel at (x, y); assumes valid fbuffer[[] array
//=========================================================================
void getRGBPIXEL(RGBPIXEL *rgb, int x, int y)
{
   int index;

   index = (linetable[y] + x) * 3;
   rgb->r = (double)fbuffer[index];
   rgb->g = (double)fbuffer[index + 1];
   rgb->b = (double)fbuffer[index + 2];
}
//=========================================================================
// Generate random paramters for the XFRM struct "xform"
// I've altered the first 3 static arrays somewhat from the original code
//=========================================================================
void random_xform(int ivar)
{
   int i, j, k, var;

   //----- number of xforms to be used in the current image -----
   static int xform_distrib[] = {2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 6, 6};

   //----- which formula(e) to use in the active xforms -----
   //----- -1 means a random mix; otherwise all xforms will -----
   //-----    use the same formula (0 - 6) -----
   static int var_distrib[] = {-1, -1, -1, -1, -1, -1, -1, 0, 1, 2,
			       3, 4, 5, 6};

   //----- weighted choices for a random mix of formulae -----
   static int mixed_var_distrib[] = {0, 0, 0, 1, 1, 2, 2, 2, 3, 3,
				     4, 4, 4, 5, 5, 6, 6};

   p.active = random_distrib(xform_distrib);
   var = (ivar < 0) ? random_distrib(var_distrib) : ivar;
   for (i = 0; i < p.active; i++)
   {
      x[i].density = 1.0 / (double)p.active;
      x[i].color = i == 0;
      for (j = 0; j < 3; j++)
	      for (k = 0; k < 2; k++)
	         x[i].c[j][k] = qsrandom11();
      for (j = 0; j < NVARS; j++)
	      x[i].var[j] = 0.0;
      if (var >= 0) x[i].var[var] = 1.0;
      else x[i].var[random_distrib(mixed_var_distrib)] = 1.0;
   }
   for (; i < NXFORMS; i++) x[i].density = 0.0;
   p.scale = 1.0;
   p.center[0] = p.center[1] = 0.0;
}
//===========================================================================
//    A pseudo-random-number generator for double-precision floats between
// 0 and 1. Can generate virtually all numbers of 9 decimal places, with
// another 9 superfluous digits tacked on.
//===========================================================================
double qsrandom01(void)
{
   int random_integer, temp_integer;
   double random_double, temp_double;

   random_integer = random(RAND_MAX);
   random_double = (double)random_integer / RAND_MAX;

/*----------------------------------------------------------------------
   At this point, regardless of the precision, random_double is limited
      to 32,767 non-consecutive possibilities. The next part of the
      routine adds a second random number between 0 and 0.000030518 to
      the first number, increasing the possibilities to 32,767 * 30,518.
      This effectively allows generation of virtually all numbers of 9
      decimal places. The next 9 places of each number are really just
      artifacts of the precision of double-precision floats.
   0.000030518 is chosen because (RAND_MAX - 1) / RAND_MAX =
   0.999969481490524004; adding 0.000030518 keeps the maximimum value
   of random_double less than 1. (0.999999999490524004)
----------------------------------------------------------------------*/

   temp_integer = random(30519);
   temp_double = (double)temp_integer / 1000000000L;
   random_double += temp_double;

   return(random_double);
}
//=========================================================================
// Return a double from -1.0 to 1.0
//=========================================================================
double qsrandom11(void)
{
   return((qsrandom01() * 2.0) - 1.0);
}
//=========================================================================
// Initialize various GUI and flame items
//    get paths to help file
//    set xres and yres
//    set parameter defaults and opening image parameters
//=========================================================================
void initialize(void)
{
   int i, j;

//   GetCurrentDirectory(MAX_PATH, helppath);
//   strcat(helppath, "\\qsflame.hlp");

   xres = yres = 600;
   oversample = 1;
   xos = xres * oversample;
   yos = yres * oversample;

   //----- parameter defaults -----
   p.iter = 100000;          // not changed for now
   p.radius = 4;             // 1 - 100
   p.white_level = 255;      // 1 - 255; not changed for now
   p.gamma = 1.0;            // 0.1 - 10.0; 1.0 == no effect
   p.brightness = 0;         // -100 - 100; 0 == no effect
   p.contrast = 1.0;         // not changed for now
   p.center[0] = 0.0;
   p.center[1] = 0.0;
   p.scale = 1.0;            // < 1.0 == zoom out; > 1.0 == zoom in
   p.pixels_per_unit = 100;  // 100; not changed for now
   p.sample_density = 5.0;   // 5.0; not changed for now

   for (j = 0; j < 2; j++)
   {
      x[j].density = 0.5;
      x[j].var[0] = 1.0;
      for (i = 1; i < 7; i++) x[j].var[i] = 0.0;
   }
   x[0].color = 1.0;
   x[1].color = 0.0;

   x[0].c[0][0] = -0.99;
   x[0].c[0][1] = -0.05861629;
   x[0].c[1][0] =  0.806310617;
   x[0].c[1][1] = -0.688782922;
   x[0].c[2][0] = -0.368144780;
   x[0].c[2][1] =  0.242713706;
   x[1].c[0][0] =  0.191625721;
   x[1].c[0][1] =  0.724478896;
   x[1].c[1][0] = -0.569627979;
   x[1].c[1][1] =  0.390423292;
   x[1].c[2][0] = -0.021881771;
   x[1].c[2][1] = -0.208493301;

   p.active = 2;
   p.scale = 0.5;
}
//=========================================================================
// Initialize more items, requiring the main window to have been created,
//    so initgraphics() will have been called and graphics mode determined
// Allocate input buffer for filter and set flag; disable filter menu
//    items if failure
// Allocate hit buffer and set flag
//=========================================================================
void initializemore(void)
{
   HMENU hMenu;
   int j;

   for (j = 0; j < yos; j++) linetable[j] = j * xos;
   setup_gamma(p.gamma);
   hMenu = GetMenu(g_hWnd);
   EnableMenuItem(hMenu, IDM_COLORPRM, MF_BYCOMMAND | MF_ENABLED);
   if (p.brightness != 0)
      EnableMenuItem(hMenu, IDM_BRIGHT, MF_BYCOMMAND | MF_ENABLED);
   if (p.gamma != 1.0)
      EnableMenuItem(hMenu, IDM_GAMMA, MF_BYCOMMAND | MF_ENABLED);

   if ((fbuffer = (unsigned char *)LocalAlloc(LPTR, xos * yos * 3)) != NULL)
   {
      fbufferflag = TRUE;
      setup_filter(p.radius);
   }
   else
   {
      oversample = 1;
      xos = xres;
      yos = yres;
      for (j = 0; j < yos; j++) linetable[j] = j * xos;
      wsprintf(string, "Error creating input image buffer.");
      MessageBox(NULL, string, "Memory Allocation Error",
                 MB_OK | MB_ICONSTOP);
      EnableMenuItem(hMenu, IDM_APPLY, MF_BYCOMMAND | MF_GRAYED);
   }
   if ((hits = (short *)LocalAlloc(LPTR, xos * yos * sizeof(short))) != NULL)
      hitflag = TRUE;
   else
   {
      wsprintf(string, "Error creating pixel hit count buffer.");
      MessageBox(NULL, string, "Memory Allocation Error",
                 MB_OK | MB_ICONSTOP);
   }
}
